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ABSTRACT 

The evolution of the luminous infrared galaxy population is explored using a pure luminosity evo- 
lution model which incorporates the locally observed luminosity-temperature distribution for IRAS 
galaxies. Pure luminosity evolution models in a fixed ACDM cosmology are fitted to submillimeter 
(submm) and infrared counts, and backgrounds. It is found that the differences between the locally 
determined bivariate model and the single variable luminosity function (LF) do not manifest them- 
selves in the observed counts, but rather are primarily apparent in the dust temperatures of sources 
in flux limited surveys. Statistically significant differences in the redshift distributions are also ob- 
served. The bivariate model is used to predict the counts, redshifts and temperature distributions of 
galaxies detectable by Spitzer. The best fitting model is compared to the high-redshift submm galaxy 
population, revealing a median redshift for the total submm population of z = 1.81q 4? m good agree- 
ment with recent spectroscopic studies of submillimeter galaxies. The temperature distribution for 
the submm galaxies is modeled to predict the radio/submm indices of the submm galaxies, revealing 
that submm galaxies exhibit a broader spread in spectral energy distributions than seen in the local 
IRAS galaxies. 
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1. INTRODUCTION 

With an energy comparable to the optical/UV back- 
ground, measurements of the far infrared background 
reveal it to peak at around 200/im (Puget et al. 1996; 
Fixsen et al. 1998), arising from dust-reprocessing of high 
energy radiation and star formation and AGN activity 
in z > galaxies. Such obscuration could be hiding ap- 
proximately half of the massive star formation activity 
over the history of the Universe (e.g. Blain et al. 1999a). 
Clearly, to unravel the cosmic history of star formation 
the evolution of the infrared galaxy population needs to 
understood. 

Studies of the evolving infrared galaxy populations 
have typically assumed a small range of template galaxy 
spectral energy distributions (SEDs), or even a single 
SED. However, infrared galaxies span a large range in 
properties, with 60 /im/100 /im colors spanning ~ 0.4 dex 
for a given luminosity (Soifer & Neugebauer 1991, Chap- 
man et al. 2003a). This distribution in dust SEDs implies 
that substantial numbers of both extremely luminous, 
yet cold galaxies, as well as low luminosity, hot galaxies 
are found. 

In a previous paper (Chapman et al. 2003a - hereafter 
C03) an evolving distribution bivariate in luminosity and 
color, $(£, C), was presented, providing consistency with 
the broad distribution of submm galaxies observed lo- 
cally. This earlier paper demonstrated that flux-limited 
surveys in various infrared and submillimeter (submm) 
bands would subsume a non-negligible fraction of both 
cold, luminous galaxies and hot, faint galaxies, the pre- 
dictions borne out in observations of low and moder- 
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ate redshift IRAS galaxies, as well as micro Jansky radio 
sources. The overall conclusion of this study was that 
surveys which select objects at either the cold Raleigh- 
Jeans tail of the dust SED, or the hot Wien tail, will pref- 
erentially detect appropriately cold or hot objects for a 
given luminosity class, in much larger numbers than ex- 
pected if the temperature distributions were not taken 
into account. 

Submm-luminous, extragalactic sources (Smail, Ivison, 
Blain 1997) currently provide our only means of study- 
ing the high redshift infrared galaxy population (Chap- 
man et al. 2003b). These galaxies are now routinely 
detected with the SCUBA/JCMT and MAMBO/IRAM 
instruments, and over 200 blank field sources are now 
cataloged from cluster lensed surveys (see Blain et al. 
2002 for a summary). The tight correlation observed 
locally between thermal FIR emission and synchrotron 
radio emission (Helou et al. 1985, Condon 1992) allows 
the identification of submm sources, recovering ~65% of 
the blank field counts (Ivison et al. 2002; Chapman et 
al. 2003a). With spectroscopic redshifts for the submm 
galaxies (SMGs), the possibility exists to assess more 
subtle effects in our model, and to understand whether 
our model can explain the range in properties subsumed 
by the SMGs, as well as explore the implication of various 
selection approaches employed in submm surveys. 

In this paper, the bivariate model is examined in light 
of the mult i- wavelength counts and backgrounds, using 
a parameterization of pure luminosity evolution. The 
model is first used to explore predictions for the Spitzer 
Space Telescope galaxy populations in terms of their red- 
shift, color, and temperature distributions. The model is 
then employed to better understand the range in prop- 
erties and redshift distributions sampled by the radio- 
identified submm galaxy population, and the impor- 
tance of the luminosity-dust temperature distribution, 



2 



CD 



m 

A 




4000 



100 



0.001 0.01 0.1 1 10 100 
S (mJy) 



1000 



X(/nm) 
1000 100 



10 



• Bivariate model 
□ Fixsen et al. (1 
1 IRAS 



□ 
□ 



10 1 



10 1 



10 1 



i/(Hz) 



Fig. 1. — The upper panel presents the bivariate model count 
from our best fit simulation at 850 jum (red), 170 /im (purple), 
160 /im (green), 70 /mi (cyan), 24 /im (blue), and 15 /mi (yellow), 
compared to SCUBA and ISO counts from the literature (points). 
For comparison, we overlay the model counts of Lagache, Dole & 
Puget (2003) as dashed lines. In the lower panel, the integrated 
bivariate model is compared to the Cosmic Infrared-Background 
of Fixen et al. (1998). The upper-limits at 12/xm and 24/ira are 
drawn from Scott et al. (2001). 



C — T, and by extrapolation the properties of the entire 
submm population. All calculations are undertaken in a 
VLq — 0.3, Ao = 0.7, h=0.65 cosmology. 

2. EVOLUTION OF THE BIVARIATE LF 

The evolutionary model is anchored to the local FIR lu- 
minosity function (LF), constructed from the 1.2 Jy sam- 
ple of IRAS galaxies (Fisher et al. 1995, C03). Recent 
work has emphasized the variation in dust temperature 
found in ULIRGs with Td as low as 25 K (Chapman et al. 
2002c). The adopted form of the LF represents the distri- 



bution in dust temperatures found in the 1.2 Jy sample, 
parametrized by the 60/mi/100/im color (a full charac- 
terization of this LF can be found in C03). 

The local bivariate LF is then evolved using pure lu- 
minosity evolution with redshift. While a range of func- 
tional forms have been employed in the literature in the 
study of the various IR, submm and radio populations, 
Blain et al. ( 1999a, b) have demonstrated that models 
which are not strongly peaked, in particular those which 
flatten beyond a certain redshift, have problems over- 
predicting the submm background. While it is clear that 
it will not be possible to distinguish minute details of 
evolutionary form, this study aims to constrain its gross 
properties. 

2.1. Fitting the multi-wavelength counts and 
backgrounds 

Studies invoking a pure density evolution of the FIR 
galaxy population have been shown to grossly over pre- 
dict the cosmic infrared background and are hence un- 
physical. Other authors have considered both density 
and luminosity evolution of the FIR population, success- 
fully recovering both the number counts and backgrounds 
(Franceschini et al. 2001; Chary & Elbaz 2001; Xu et al. 
2004). However, it has been shown that these observa- 
tional properties can be explained within models adopt- 
ing only luminosity evolution (Blain 1999a). Given its 
simplicity (with reduced parameters), we adopt a similar 
pure luminosity evolution to examine the cosmic history 
of FIR galaxies. 

The evolution is modeled in a simple form; $(L,C) = 
$>o(L/g(z), C). The color term of our LF (C f rest frame 
60/im/100/im flux ratio) does not evolve, but rather 
continues to scale with FIR luminosity as found locally. 
C03 and Blain, Barnard & Chapman (2003) have demon- 
strated that this lack of significant 60/im/100/im color 
evolution appears to hold at least out to moderate red- 
shifts (0.3 < z < 0.9) in a sample of IRAS detected 
ULIRGs of Stanford et al. (2000). As the form of the 
C — T distribution remains fixed for all redshift s, the 
only remaining parameter is g(z). It should be noted that 
Dunne et al. (2000) suggested that the IRAS galaxy lu- 
minosity function may be incomplete, missing cold galax- 
ies. Recently, this conclusion has been strengthened by 
observations by Klaas et al. (2001) and Bendo et al. 
(2002; 2004). Given that the current study is tied to 
IRAS, this effect will accentuate the features of cold 
galaxies presented in this paper, although the degree of 
IRAS underrepresentation of cold galaxies must be quan- 
tified to fully explore this. 

Adopting a Monte Carlo approach, luminous infrared 
galaxies are drawn randomly from the evolving distribu- 
tion function. A galaxy thus selected from this model 
Universe is then assigned a template spectral energy dis- 
tribution from the catalog of Dale et al. (2001, 2002), 
parameterized by the 60/im/100/im color, and normal- 
ized to the FIR luminosity. This model scenario does not 
incorporate the intrinsic scatter in the FIR/radio corre- 
lation (0.2 dex), as was done in Chapman et al. (2002b; 
hereafter C02). Rather this study concentrates on the 
properties of the intrinsic dust temperature distribution 
which locally show a larger scatter (0.3 dex) than the 
FIR/radio correlation, and should be expected to domi- 
nate the high redshift radio and submm properties. 
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Fig. 2. — Redshift distributions at 24 /im, 70 /im, 160 jum, and 850 jum. The distributions are shown for both a deep (left) and a shallow 
(right) survey, in the bivariate (solid line) and single variable (dashed line) models. The effect of the bivariate distribution is to increase 
the numbers of both high and low-redshift sources relative to the single variable model, due to the temperature spread. The shallow survey 
is comparable to the expected FLS survey, with the 850 jum survey being typical of current wide-field SCUBA surveys. The deep survey is 
scaled at all wavelengths by a factor of 5. 
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Fig. 3. — Color-color diagrams of galaxies in a 24 /im flux limited survey (S24yum > 0.2 mJy). The left-hand plot presents 70/ira/160/ira 
versus 8/ira/24/i,m for the univariate case, whereas the right-hand panel presents the same distribution for the bivariate case. Objects are 
color-coded by redshift: < z < 1 (cyan), 1 < z < 2 (blue), 2 < z < 3 (green), z > 3 (red). 



C02 adopted an evolution function with a simple power 
law peak, g(z) oc (1 + z) A out to a break redshift, 
and dropping thereafter as g(z) oc (1 + z) -4 , with no 
discontinuity in the zeroth derivative. This power-law 
index was chosen provided reasonable descriptions of 
submm and radio data, allowing a fit to the sub-mm 
counts and background (C02). This evolution form 
coupled with the bivariate model provided a reason- 
able fit over the submm fluxes represented by our data 
(5<S85o,um<15mJy). However, the model produced too 
many bright submm sources, with an increasing excess 



of >20mJy sources with increasing peak redshift. This 
toy model cannot be considered physical, and the culprit 
for the excess bright sources is found to be in the peak 
turnover in g{z) producing a spike of luminous sources 
at the tip of our color-magnitude diagram (CMD) cou- 
pled with the range of Td for each luminosity. While 
this peak model was useful for illustrating various se- 
lection effects (as in C03), it must be replaced with a 
more realistic evolution form to provide a more accu- 
rate description of the submm population and, therefore, 
the functional form of Blain et al. (1999b) is adopted 
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g(z) = (1 + z) 3/2 sech 2 [b\n(l + z) - c]cosh 2 c. 

In constraining the form of the above evolution model, 
galaxy submm/radio distributions with varying values of 
b and c were constructed and compared to the differen- 
tial number count models at 850/xm (Barger, Cowie & 
Sanders 1999) and 15/xm (Elbaz et al. 1999). For partic- 
ular set of values of b and c, the number counts are scaled 
to minimize the residuals between the Monte Carlo simu- 
lation and the modeled number counts. For the purposes 
of this study, it was assumed that the models possessed 
an uncertainty of 0.2 dex and the resulting best values 
of b and c are 2.10±g"Jg and 1.81±g;|| respectively. The 
range of parameters effectively shift the %> ea k of the peak 
evolution function such that z pea k = 1.8^4. 

Bivariate model count from our best fit simulation at 
850 //m, 170 //m, 160 /xm, 70 /xm, 24 /xm, and 15 /xm, are 
shown in the upper panel of Fig.^ The models are over- 
laid on measured SCUBA (850 /xm- Blain et al. 2002) 
and ISO (15 /xm - Elbaz et al. 1999; and 170 /xm- Dole 
et al. 2001) counts from the literature. For compari- 
son, we overlay the model counts of Dole et al. (2003) 
as dashed lines, which agree remarkably well with our 
models. At the brightest end of the counts, our model 
is constrained by the length of our Monte Carlo runs, 
and we are dominated by small number statistics. In the 
lower panel of Fig. [l] presents compares the integrated 
flux in the bivariate model, compared to the observed 
Cosmic Infrared Background (CIB; Fixen et al. 1998); 
again, the best fit model from the counts accounts for 
the observed CIB distribution. 

Previous modeling (C02) found a reasonable fit to the 
data with a peak evolution redshift of z = 3, when 
all objects were assumed to have hot dust temperatures 
(~ 50 K), similar to the local ULIRG, Arp220. The cur- 
rent model ties the dust temperature of a source to its 
luminosity as observed locally for the IRAS 1.2 Jy sam- 
ple galaxies, whereby a galaxy with a FIR luminosity 
of 10 12 L has Td=35K and dust emissivity /3 = 1.6, 
as described in Dale & Helou (2002). In addition, 
the 60/xm/100/xm color distribution observed locally is 
adopted in the FIR luminosity function, resulting in a 
tail of cold luminous galaxies which are preferentially se- 
lected with SCUBA at 850 /xm (e.g., Blain 1999b, Eales 
et al. 1999, Chapman et al. 2002c). As expected (see 
also the discussion in C02), the best-fitting peak evolu- 
tion redshift must be lower if a population of luminous 
colder sources exists at high redshifts. 

3. PREDICTIONS FOR SPITZER AND SUBMILLIMETER 
GALAXIES 

3.1. Redshift & color- color distributions 

Taking our best-fit evolution functions for both the 
bivariate and single variable models, we predict redshift 
distributions for galaxy surveys at 24 /xm, 70 /xm, 160 /xm, 
and 850 /xm, covering an area of one square degree. We 
place flux density thresholds on our simulation, compa- 
rable with a shallow Spitzer observation such as the First 
Look Survey (http://ssc.spitzer.caltech.edu/fls). 
A deep survey is also analyzed, representing flux densi- 
ties five times fainter in all bands (ignoring the effects 
of confusion in the longer wavelength beams for now). 
For the submm surveys, the shallow survey is set at the 
confusion limit of SCUBA (2mJy). The survey reach- 
ing five times this depth will be achievable with future 



instruments such as ALMA. 

The comparisons are shown in Fig. [2j Due to the 
spread in temperature at each luminosity, with increas- 
ing numbers of hot sources lying at higher redshifts and 
conversely colder sources at lower redshift, the effect of 
the bivariate distribution is to increase the numbers of 
both high and low-redshift sources relative to the sin- 
gle variable model. In considering photometric redshifts 
for Spitzer sources, or indeed simply trying to pre-select 
sub-populations of Spitzer sources with a given range in 
redshift or temperature properties, it is, therefore, crucial 
to have a calibrated estimate of the source distributions 
in the various color-color planes. 

Color-color diagrams of galaxies in a 24 /xm flux lim- 
ited survey (S24 Mm > 0.2 mJy) are shown in Fig. 03 The 
left-hand panel presents the distribution of galaxies as- 
suming the univariate luminosity-temperature relation, 
whereas the right-hand panel presents the correspond- 
ing distribution for the bivariate case. In each panel, 
objects are color-coded by redshift: < z < 1 (cyan), 
1 < z < 2 (blue), 2 < z < 3 (green), z > 3 (red). It 
is apparent that the distribution of galaxies in the two 
figures is markedly different, with the introduction of a 
bivariate form of the luminosity function changing the 
relative density of galaxies in the color-color plane; this 
behavior is also seen in other color-color planes. Clearly, 
while different redshift regimes can be delineated in the 
color-color plane, allowing the determination of photo- 
metric redshifts for galaxies, the redistribution of galax- 
ies in the plane due to the introduction of the bivariate 
distribution implies that such determinations depend on 
the nature of the luminosity-temperature relationship. 

While our model has not yet been tuned to fit the deep, 
mult i- wavelength Spitzer counts, it already emphasizes 
that the source distributions as a function of redshift and 
SED can vary significantly, from the bivariate case to the 
single temperature luminosity models. Upcoming papers 
will explore this issue more fully, once deep Spitzer data 
is available to further constrain our models. 

In C03 we presented the baseline assumption in this 
model: the local Ssso^m/Si^GHz as a function of C holds 
at all redshifts, simply extrapolating the local C — C 
relation to the required luminosities. We have now fit 
this model explicitly to the pie- Spitzer submm and mid- 
infrared counts to constrain the evolution function, and 
used this tuned model to make predictions on the popu- 
lations observed in the Spitzer wavelengths (8 /xm, 24 /xm, 
70 /xm, and 160 /xm). 

In Fig. EI we show the Ltir distribution as a function 
of redshift. Each pair of panels represents a particu- 
lar waveband (8 /xm, 24 /xm, 70 /xm, and 160 /xm respec- 
tively), with the left-hand panel in a pair representing the 
bivariate model (blue markers) and the right-hand panel 
presenting the univariate case (red markers). Until we 
apply a flux limit to the figure (darker points), there is 
no difference between the visualizations since they have 
the same luminosity evolution formalism. Fig. 0] is our 
Monte Carlo representation of the evolving LF; vertical 
slices reveal the dual power law $(£, C) at each redshift. 

When flux limited surveys are considered in the context 
of Fig. 01 differences in the two models become manifest. 
In the case of the single variable distribution, each Ltir 
point maps uniquely to a flux for a given wavelength. 
However in the bivariate LF, each Ltir point corresponds 
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Fig. 4. — Comparison of our Monte Carlo representation of the evolving C) (blue points) to an evolving univariate model LF (red 

points) with a one-to-one mapping of luminosity to color. Each point represents only luminosity and thus the distributions are equal in 
both cases. The difference between models is apparent when flux limited surveys are defined on the model points. We overlay flux limited 
samplings of the model points with darker symbols. In the single variable LF model, luminosities map uniquely to fluxes and lines can be 
drawn which reflect fixed temperatures. The upper panels are for 8 /xm and 24 /xm, while the lower panels are for 70 /im and 160 /im. 
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Fig. 5. — Comparison of the dust temperature distributions as a function of redshift for the <E>(£, C) (solid lines) to an evolving model 
LF with a one-to-one mapping of luminosity to color (dashed lines). Distributions are shown for hot (red lines - Td>40K), medium (green 
lines - 28<Td<33K), and cold (blue lines - Td<25K). The panels are divided by wavelength, as indicated. 



to a probability distribution of fluxes corresponding to 
the log-normal distribution in S850/xm/Si.4GHz and the 
associated range of SED templates that can be tied to 
the Ltir value. The sensitivity limits in the various 
bands are shown for the deepest surveys with Spitzer. 
The structure in the 8 /im and 24 /im surveys is a result 
of PAH bands (rest ~10 /im) being redshifted through 
the Spitzer 24 /im filter. At 70 /im and 160 /im, the SED 
is relatively smooth. 

The effect is subtle in Fig. as both the luminos- 
ity function and the color distribution are scattering the 
observed fluxes, largely canceling dramatic differences in 
the effective luminosity limits probed with redshift. How- 
ever, the most important difference between the bivariate 
and luminosity-only models is apparent in Fig. 0] in the 
simpler $(£) model, the flux limit for a given wavelength 
translates at each redshift into a transition range of lu- 
minosities, within which galaxies are or are not detected 
depending on their color. For surveys selecting sources 
along the hot dust side of the grey-body peak, that be- 
ing the case for all the accessible wavelengths of Spitzer 
except 160 /im, colder luminous sources will be missed 
and hotter low luminosity sources will be preferentially 
detected. 

The sensitivity limit of Spitzer is adversely affected by 



the steep Wien slope of the far-mid-infrared SED, mak- 
ing more distant sources difficult to detect and resulting 
in the steep sensitivity curve in Fig. 21 The result is that 
near the sensitivity limit of Spitzer, surveys will be dom- 
inated by the large numbers of sources lying both above 
and below any fixed temperature cut. For hotter dust 
temperatures, an excess of a factor greater than two of 
sources is predicted by the C) model. All sources 

are luminous enough to be detected regardless of their 
temperature, and there are far more lower luminosity 
sources which are boosted by the bivariate LF to hotter 
dust temperatures than vice versa. 

Fig. demonstrates these model differences explicitly, 
demonstrating that many more sources are detected in 
the extremes of the temperature distributions in the 
C) model. Moreover, the predicted distributions are often 
radically different in form in the C) case. The dou- 
ble peaked profile apparent in the redshift distributions 
presented in Figure [21 is a result of temperature differ- 
ences in the underlying IR population. Fig.0 therefore, 
illustrates the strongest differences between the univari- 
ate and bivariate models. 

Recently deep counts at 24/im, 70/im and 160/im, ob- 
tained with the MIPS intrument on Spitzer, have be- 
come available (e.g. Dole et al. 2004; Papovich et al. 
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Fig. 6. — As for the upper panel in FigurefTl but now the bivariate 
model is compared to the recently compiled Spitzer counts at 24/ira 
(blue), 70fim (cyan) & 160/xra (green) (Dole et al. 2004; Papovich 
et al. 2004). 



500 



a 100 




redshift 

Fig. 7. — Ss50)Ltm/Si.4GHz indices for SMGs with redshifts, com- 
pared to our model prediction, with the dashed lines corresponding 
to lcr scatter in the population of model galaxies. Note that the 
model includes only the range of dust temperatures observed lo- 
cally in the IRAS population, whereas a variation in the radio to 
farlR ratio likely increases the scatter further. 



2004). This presents further opportunities to test the 
efficacy of the bivariate model. Figure |S| presents the 
bivariate model prediction for the Spitzer MIPS bands 
with the observed counts overlaid. Generally, the ob- 
served trends are reproduced, but discrepancies are ap- 
parent, most notably a deficit in the Spitzer counts be- 



tween lmJy - 10m Jy. Addressing these discrepancies is 
beyond the scope of this present work and will be re- 
served for further study. 

4. STUDYING THE SUBMILLIMETER GALAXIES 

Spitzer data is just beginning to be transmitted back to 
Earth, and analysis and spectroscopic followup will take 
a significant effort. Detailed comparison of our models 
with the Spitzer sources will be the focus of an upcom- 
ing paper, once an initial census of the Spitzer surveys 
are complete. However our understanding of the deepest 
850 fim submm sources has recently reached a relatively 
mature state, with the measurement of spectroscopic red- 
shifts for a large sample (Chapman et al. 2003b; Chap- 
man et al. in preparation). In this section, we compare 
our predictions for the submm galaxies (SMGs) directly 
with the measurements. 

The S850/xm/S 1.4GHz ratio provides a measure of the 
degenerate quantity (l+z)/Td (Blain 1999b), coupled 
with any evolution in the Far-IR/radio correlation. As 
the redshift parameter has now been independently 
measured for the SMGs, we can directly compare the 
S850/im/Si.4GHz predictions of our model to the SMGs. 
Fig. shows the S850/zm/Si.4GHz indices for SMGs with 
redshifts, compared to our model prediction. Note that 
the model includes only the range of dust temperatures 
observed locally in the IRAS population, whereas a varia- 
tion in the radio to farlR ratio likely increases the scatter 
further. 

The actual submm population appears to span a larger 
range than the local IRAS galaxy population, although 
there are some caveats to consider. The SMG surveys 
are missing cold sources in the radio detection criterion, 
and thus there is some asymmetry in the actual observa- 
tions towards hotter dust temperatures. At low redshifts, 
there are some apparently very cold galaxies, which are 
not predicted by the local IRAS extrapolation. If the 
SMG identifications are correct, they may signal a rare 
population which has no local analogs. 

A remaining question is the interpretation of sources in 
our model. Ongoing debates discuss whether accretion 
power from super-massive black holes or bursts of star 
formation is heating the dust in ULIRGs. Our model 
has not assumed either explicitly. Rather, we have taken 
the distribution of 60/im/100/im color found locally in 
the IRAS galaxy population and assumed it describes 
galaxies at all redshifts. However, AGN are known to 
typically have warmer colors, and our model may not 
accurately reflect their possibly larger contributions at 
higher redshift, higher luminosity. 

A test for AGN contribution is whether the SMGs are 
unusually bright in the X-ray, where even obscured AGN 
should reveal themselves, penetrating high columns of 
gas. The radio identification of the submillimeter sources 
has allowed their location to be pinpointed. In Barger 
et al. (2001a, b), it was pointed out that approximately 
10% of the X-ray sources were submillimeter detected 
in the Chandra Msec exposure centered on the HDF. 
The 2 Msec exposure in this field has detected the major- 
ity (~ 90%) of the radio-SMGs (Alexander et al. 2003). 
As this radio identified submillimeter population repre- 
sents ~65% of the total blank-field SCUBA population 
at S , 850/im>5 mJy, this implies a fraction > 60% of the 
total SCUBA population is detectable in a 2 Msec Chan- 
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Fig. 8. — Predicted redshift distributions for radio (dashed) and submm (solid) wavelengths. The left panel shows the prediction for the 
bivariate model, while the right panel shows the single variable LF model. Mild differences in the predicted N(z) are observed between 
models, which the data cannot differentiate between. The dark shaded region in both panels reveals the submm sources which are not 
expected to be detectable in the radio. The light shading shows the spectroscopic desert where no strong emission lines are shifted into 
the optical atmospheric window. Overlaid as a histogram is the measured redshift distribution for radio detected submm galaxies from 
Chapman et al. (2003c). 



dra exposure. The typical inferred star formation rate 
from this analysis is roughly 1000 M yr _1 , similar to 
that deduced from the submillimeter luminosities. 

The sources with considerably higher X-ray luminosi- 
ties, where the X-ray implied SFR would severely ex- 
ceed the submillimeter estimate, suggest an AGN must 
be generating the X-ray emission, and we emphasize that 
our model does not account for such sources explicitly. 
Note however that in Fig. only one source explicitly 
stands out from the distribution due to radio excess over 
the model expectations. 

Let us now look specifically at the differences between 
the bivariate and single variable models. The counts 
and redshift distributions are natural measurable outputs 
from our model. The univariate model predicts a one-to- 
one mapping of S850/zm/Si.4GHz with redshift [for a fixed 
Q-value (Helou et al. 1985)]; clearly the observed scatter 
seen in Figure [3 cannot be explained such a simple map- 
ping. The bivariate model naturally introduces a scatter 
into this relationship, and the dashed lines in Figure [3 
correspond to the la width of the Ssso^m/Si^GHz distri- 
bution as a function of redshift. It is a straw man argu- 
ment to compare a single SED model to our data, which 
provides an unphysically narrow range of submm/radio 
colors (e.g., C02). A more realistic approach is to con- 
strain the width of the color distribution of our best fit- 
ting bivariate model until there is only a single tempera- 
ture SED associated with a given luminosity. This single 
variable model can be thought of as mapping the dust 
temperature monotonically to the source luminosity. 

We first extract the counts from our models and com- 
pare directly with the measured field counts. Both mod- 
els are able to fit the counts well within current mea- 
surement errors. The submm counts are insensitive to 
the width in color, C , in our model due to the degen- 
eracy in Td versus redshift. Therefore neither of these 
models is better constrained by the submm counts. The 
far-IR background (Puget et al. 1996; Fixsen et al. 1998) 
also provides little constraint on the detailed form of the 
high-z evolution, requiring only that the evolution func- 



tion fall off quickly enough at high redshifts to avoid 
generating too high an energy density. Both models 
generate submm and far-IR background values midway 
between the measurements of Puget and Fixsen. The 
model is therefore self-consistent as a generalization to 
the evolution of the entire submm population. While 
the model does not provide a unique description of the 
submm galaxies, it is physically motivated by the de- 
tailed properties of local and moderate redshift IRAS 
galaxies. 

We can use our model to predict the redshift distri- 
butions for radio (1.4 GHz) and submm (850 /im) galax- 
ies. In Fig. [HI we show the prediction for the bivariate 
model (left panel) compared to the single variable LF 
model (right panel). The bivariate scenario marks a dif- 
ference in the predicted N(z), where a deficit of sources 
at the redshift peak are seen relative to the single vari- 
able model, with a corresponding boost to the higher 
redshift sources. This behavior was described in C02, 
where the far-IR/radio distribution provides a similar ef- 
fect. The shaded regions between the radio and submm 
model curves reveal the submm sources which are not 
expected to be detectable in the radio. In the bivari- 
ate model, a range of hotter dust temperatures for a 
given luminosity lead to a higher redshift range for ra- 
dio detect ability. While not statistically significant, the 
bivariate models appears to follow the observed redshift 
distribution. 

Understanding of the full redshift distribution for the 
submm population through deep optical and millimeter 
spectroscopy will provide a strong test of these models 
and the applicability of the local bivariate C — T relation 
to this high redshift population. The discovery of the red- 
shift distribution for the radio detected submm galaxies 
(Chapman et al. 2003a, c) is shown overlaid on our mod- 
els in Fig. [HI A remarkable agreement is found between 
the data and our model, with a preference for the flat- 
tened distribution and expanded range of redshifts in the 
bivariate picture. We conclude that submm galaxies ex- 
hibit at least as large a range in dust temperatures as 
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local IRAS galaxies, and likely a larger range. 
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